function  [ar,av,phr]=fun_gdmot(RD,IPS,VSVP,AIN) 
%{
C+ SUBROUTINE GDMOT(RD,IPS,VSVP,AIN,AR,PHR,AV)
C
C  GDMOT GIVES GROUND DISPLACEMENT FOR INCIDENT P (IPS=1) OR S (IPS=2)
C  BULLEN (1963) P. 129 FOR INCIDENT P (BUT AIN=90-E)
C  AIN IS EMERGENCE ANGLE IN DEGREES RELATIVE TO VERTICAL
C  VSVP=VS/VP AT THE SURFACE
C  AV AND AR ARE THE VERTICAL AND RADIAL GROUND DISPLACEMENTS
C  FOR SV INC. BEYOND CRITICAL, PHV=PHR+90.
C     (AR MAY BE NEGATIVE IN THIS CASE - NOT REALLY AN AMPLITUDE)
C-
%}
      
      phr=0.0;
     if (AIN ==0.0) 
        ar=2.*(IPS-1);
	av = 2*(2 - IPS);
     else
        if  (IPS==1) 
          A = AIN/RD;
          B=asin(sin(A)*VSVP);
          COTANB = 1.0/tan(B);
      COTCOTB=1.-COTANB^2;
          DEN=2.*cos(A)/(sin(B)^2*(4.*cot(A)*COTANB+COTCOTB^2));
          ar=2.*COTANB*DEN;
          av=-COTCOTB*DEN;
          phr=0;
        end
       if (IPS== 2) 
          B=AIN/RD;
           COTANB = 1.0/tan(B);
      COTCOTB=1.-COTANB^2;
          SA=sin(B)/VSVP;
         if(SA <=1.0) 
            A = asin(SA);
            DEN=2.*COTANB/(sin(B)*(4.*cot(A)*COTANB +COTCOTB^2));
           av = -2*cot(A)*DEN;
            ar=-COTCOTB*DEN;
         else
            CTA=sqrt(1.-1./(SA*SA));
            DEN=2.*COTANB/(sin(B)*sqrt((4.*CTA*COTANB)^2 +COTCOTB^4));
            av=-2.*CTA*DEN;
            ar=0.;
            phr = 90.0;
               if (COTCOTB>0.0 || COTCOTB<0.0 ) 
                ar=-COTCOTB*DEN;
               phr=atan2(4.*COTANB*CTA,-COTCOTB^2)*RD;
               end
         end
        end
     end
     
end


